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1.  Introduction 

The  Quantitative  Feedback  Theory  (QFT)  technique  developed  by  Isaac  Horowitz  over  a 
number  of  years,  is  perhaps  the  only  controller  design  methodology  that  enables  a  controller  to 
be  designed  to  a  given  specification  in  a  transparent  quantitative  manner.  By  this  is  meant  that 
there  is  a  definite  quantitative  measure  of  the  closeness  of  the  design  to  an  optimum.  A  major 
advantage  of  QFT  is  the  fact  that  the  trade-offs  between  the  constraints  and  the  set  of  design 
criteria  are  visible  to  the  designer  in  a  transparent  manner  at  all  stages  during  the  actual  design  I. 
process,  rather  than  at  the  end,  as  is  the  case  with  ‘black  box’  synthesis  techniques  such  as  hW 
or  LQG  optimal  control.  ^ 

The  manual  QFT  method  introduced  by  Horowitz  and  others  in  1972''"f3T''6f  represented 
a  major  breakthrough  in  the  quantitative  design  of  robust  controllers.  However,  the  method 
is  extremely  labour  intensive  and  the  final  loop-shaping  stage  of  the  design  process  requires 
substantial  practice  and  expertise  and  it  is  believed  that  for  this  reason,  the  method  has  not 
been  as  widely  accepted  as  it  deserves  to  be. 

This  report  details  research  carried  out  to  develop  a  computer-based  method  for  optimal 
loog-shaping  in  QFT.  Although  some  work  has  already  been  done  in  this  area  by  Gera  and 
Horowitz''[ij>|n  1980,  no  practical  implementation  details  have  been  published.  We  believe  that 
in  OptComp  we  have  made  good  progress  in  developing  a  program  that  enables  the  engineer  to 
use  QFT  methods  to  design  a  compensator  (or  controller)  iteratively  to  any  d^red  order,  while 
remaining  transparent  at  all  times  about  what  trade-offs  are  necessary.  — , — . 

At  each  stage,  the  designer  chooses  the  desired  form  of  compensator,  and  the  points  in 
the  Nichols  chart  through  which  it  should  pass.  The  program  fits  the  best  curve,  and  clearly 
displays  the  difference  between  the  accuracy  required  and  the  accuracy  attained,  both  graphically 
and  in  tabular  form.  The  designer  can  then  use  his  or  her  own  judgement  to  decide  whether 
a  satisfactory  result  has  been  reached  for  a  compensator  of  this  order,  and,  if  so,  whether 
a  compensator  of  higher  order  should  be  attempted.  Thus,  while  all  the  drudgery  has  been 
removed  from  the  design  process,  the  designer  still  has  the  satisfaction  of  being  closely  involved 
in  the  art  of  design  (for  example,  finding  troughs  in  the  boundaries)  and  in  assessing  trade¬ 
offs  between  accuracy  at  different  w-values,  or  between  overall  accuracy  and  the  order  of  the 
controller.  This  means  that  the  benefits  of  QFT  axg  retained,  and  its  drawbacks  removed. 


2.  Objectives  and  contractual  details 


In  terms  of  the  original  letter  of  5  May  1996,  the  contract  was  in  two  parts: 

a.  To  develop  computer  algorithms  to  enable  computer-aided  design  of  the  optimum  loop 
transmission  for  QFT  design  of  multivariable  control  systems. 

b.  To  incorporate  the  optimal  loop  transmission  algorithms  into  the  MIMOQFT  design  package 
developed  at  the  US  Air  Force  Institute  of  Technology. 

In  the  formal  contract  the  wording  is  somewhat  different:  to  undertake  preliminary  development 
of  algorithms  for  optimum  computer-aided  loop-shaping  in  control  system  design  using  QFT, 
with  no  mention  of  incorporation  into  MIMOQFT. 

The  starting  point,  as  defined  by  Professor  Houpis  in  a  meeting  at  AFIT  in  October,  is 
with  a  given  nominal  plant  Po{ju}),  a  given  set  of  w-values  wi, . . ., and  a  set  of  composite 
boundaries  BoQwi), . . .,  Bo{jijJn)  [4,  p.  12].  The  novelty  and  originality  of  OptComp  is  is  that 
it  allows  the  user  to  design  the  compensator  G{j(jj)  directly,  rather  than  obtaining  it  indirectly 
as  the  quotient  of  LoQ'w)  and  Po{j^)  [4,  p.  14]. 

The  preliminary  programs  are  written  in  Matlab,  because  of  its  power  and  simplicity  (any 
deficiencies  in  precision  are  unimportant  at  the  development  stage),  and  because  neither  of  us 
has  convenient  access  to  Mathematica  at  the  moment.  Once  the  programs  have  been  final¬ 
ized,  translation  into  Mathematica  and  incorporation  into  the  MIMOQFT  package  should  be 
relatively  minor  tasks. 


3.  Mathematical  description 

It  is  assumed  that  the  reader  is  familiar  with  at  least  the  basic  principles  of  QFT,  as  given, 
for  example,  in  [2]  and  [4].  The  conventional  process  of  compensator  design  (or  loop-shaping)  is 
clearly  described  in  [4,  p.  14],  and  consists  in  finding  the  function  G{ju;)  so  that 

W 

Lo{ju)  =  PQ{j(J)G{juj),  where  G{ju>)  =  KkGk{joj) 

fc=0 


[4,  p.l4],  equation  1.12.  However,  this  process  is  indirect,  in  that  the  designer  actually  constructs 
Lo{s),  and  “once  a  satisfactory  Lo{s)  is  achieved  then  the  compensator  is  given  by  G{s)  = 
Lois)/Po{sy  [4,  p.l4]. 

We  here  employ  an  original  approach  of  determining  the  compensator  G{s)  directly,  one 
factor  KkGk{s)  at  a  time,  until  the  desired  order  has  been  reached  or  the  compensated  function 
meets  the  required  specifications.  The  method  depends  on  two  principles: 

(a)  The  Nichols  chart  is  essentially  logarithmic  (see  Appendix  1),  so  nic(Lo(jw))  =  nic(Fo(j<^))+ 
nic(G(ia;)).  This  means  that  the  the  requirement  that  Lo{ju})  lie  on  the  boundary  Bo{ju))  is 
equivalent  to  requiring  nic(G(ja;))  to  lie  on  the  shifted  boundary  mc{Bo{joj))  —  nic(Po(jA;)) 
in  the  Nichols  chart. 

(b)  G{s)  must  be  a  rational  function  with  pole  excess  greater  than  or  equal  to  zero.  By  the 
Fundamental  Theorem  of  Algebra,  every  polynomial  with  real  coefficients  can  be  factorized 
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as  a  product  of  linear  and  quadratic  factors.  It  follows  that  G(s)  can  be  written  as  a  product 
of  rational  functions  of  the  six  standard  forms  listed  in  the  table  below,  and  therefore 

W 

nic(G(s))  =  ^nic(A'fcGjt(s)), 

*:=0 


where  each  KkGk{s)  is  of  one  of  the  six  standard  forms. 


1. 

2. 

3. 

4. 

5. 

6. 


constant/linear 
linear/linear 
linear/quadratic 
quadratic/quadratic 
constant  over  quadratic 


K/{s+b) 

K{s  +  a)/{s  +  b) 

K{s  +  a)/{s'^  +  bis  +  bo) 
K{s‘^  +  ais  +  ao)/(s^  +  i»is 
K/ +  6is  +  bo) 

K/s 


(1st  order  lag) 

(1st  order  lead-lag) 

6o)  (2nd  order  lead-lag) 
(2nd  order  lag) 
(pure  integrator). 


The  six  standard  forms  of  compensator. 


From  the  two  observations  (a)  and  (b)  above  comes  the  fundamental  result,  on  which  the  program 
is  based. 

Theorem.  Any  compensator  can  be  achieved  by  iteratively  finding  standard  compensators 
KkGk{s)  for  k  =  1,2,  ...,u;  such  that  mc(KkGk(ju)))  lies  as  close  as  possible  to  the  shifted 
boundary  mc{Bo{jijj))  —  mc{Pk-i{jo:))  for  all  given  values  of  u).  Here  Pk-i{s),  for  k  >  1, 
denotes  the  partially  compensated  plant  Fo(s)  Hr^i  f^rGr(s). 

In  practice,  the  values  wi , . . . ,  a>„  are  given,  as  well  as  the  plant  values  Pq  (j‘^i )  i  •  •  •  i  (J‘^n) , 
and  the  composite  boundaries  Bo(joJi), . . . ,  Bo(j<^n)-  The  first  step  is  to  find  a  standard  com¬ 
pensator  i^iGi(s)  so  that  nic(/fiGi(jWr))  lies  as  close  as  possible  to  the  shifted  boundary 
nic(Bo(i^r))  -  nic(Po(iWr))  for  r  =  1, . . . ,  n.  Qualitative  knowledge  of  the  shape  of  the  Nyquist 
plots  of  the  standard  compensators  in  the  Nichols  chart  is  therefore  essential,  and  representative 
samples  are  illustrated  in  Figure  3.1. 

For  minimum  phase  compensators,  the  plots  all  start  with  a  phase  angle  of  0°  when  a;  =  0, 
and  curve  anticlockwise  in  the  Nichols  chart  as  oj  increases.  For  the  lead-lag  compensators  with 
zero  pole  excess  (forms  2  and  4),  the  curves  end  with  a  phase  angle  of  0°  as  well,  the  difference 
being  that  form  2  curves  have  phase  of  one  sign  only,  whereas  form  4  curves  can  have  both 
positive  and  negative  phase,  and  cross  the  zero  phase  line  at  an  intermediate  value  of  w.  The 
curves  for  compensators  with  pole  excess  of  one  (forms  1  and  3)  all  approach  the  —90°  phase 
line  as  w  — >■  oo.  They  differ  only  in  that  the  phase  decreases  monotonically  from  0°  to  —90°  in 
form  1  curves,  whereas  in  form  3  curves  the  phase  first  increases,  then  decreases  beyond  —90°, 
and  increases  again.  Finally,  form  5  curves  approach  the  —180°  phase  line  as  a;  — >  oo,  and 
form  6  plots  lie  on  the  —90°  phase  line,  and  do  not  require  illustration. 
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Figure  3.1.  Nyquist  plots  of  standard  minimum  phase  compensators:  (a)  Form  1, 
(b)-(c)  Form  2,  (d)  Form  3,  (e)-(h)  Form  4,  (i)  Form  5.  The  circle  corresponds  to 
=  0  and  the  star  to  lj  oo. 


4.  User’s  guide  to  OptComp 

The  program  runs  under  Matlab,  even  under  the  student  edition,  if  the  boundary  matrix  is 
not  too  large.  Once  in  Matlab,  before  the  designer  can  use  OptComp,  the  following  data  must 
be  available: 

•  omega  —  a  row  vector  of  n  values  of  w,  preferably  in  ascending  order. 

•  pleuit  —  a  row  vector  of  n  nominal  plant  transfer  function  values  (in  the  Nichols  chart)  at 
s  =  ju},  one  for  each  value  of  w  stored  in  omega. 

•  bndry  —  an  mX  n  matrix  of  composite  boundary  points  in  the  Nichols  chart,  made  up  of  m 
boundary  points  for  each  of  the  n  values  of  u.  Note  that,  for  convenience,  each  boundary  has 
the  same  number  of  points;  this  can  be  achieved,  if  necessary,  by  interpolation,  repetition,  or 
omission. 

Three  sample  data  sets  have  been  provided,  and  any  one  can  be  used  (by  entering  OptComp  1, 
0ptComp2,  or  OptCompS)  for  immediate  trial  of  the  main  program.  All  important  output  is 
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written  into  a  diary  file,  which  will  be  the  default,  named  simply  diary,  unless  another  file 
name  is  specified  with  the  diary  command. 

After  entering  OptComp  at  the  Matlab  prompt,  the  user  is  reminded  of  the  facts  in  the 
preceding  paragraphs,  and  given  an  option  to  exit.  If  he  or  she  continues,  the  date  and  starting 
time  are  displayed,  followed  by  a  table  showing  the  minimum  compensation  required  (i.e.,  the 
magnitude  and  phase  of  the  shortest  distance  from  the  nominal  plant  to  the  boundary)  at  each 
w-value.  The  information  is  then  repeated  graphically  in  the  Nichols  chart  in  the  conventional 
way,  giving  the  Nyquist  plot  of  the  plant  (in  blue),  with  the  points  at  the  given  w- values  marked 
by  plus  signs,  and  the  boundaries  drawn  in  white.  The  boundaries  and  plant  points  are  each 
labelled  with  the  corresponding  w-value.  In  order  to  keep  the  boundaries  continuous,  the  Nichols 
chart  is  extended  horizontally,  if  necessary,  beyond  the  ±180°  phase  angles.  (The  user  can  press 
any  key  to  continue  the  program  whenever  it  pauses  for  inspection  of  graphical  or  other  display.) 

As  mentioned  above,  the  originality  of  OptComp  is  that  it  enables  the  designer  to  concen¬ 
trate  directly  on  the  compensation  required  at  each  w-value.  The  following  displays  all  concern 
the  Nyquist  plot  (in  the  Nichols  chart)  of  the  compensator  that  is  to  be  found.  The  first  diagram 
shows  the  differences  between  the  composite  boundaries  and  the  corresponding  plant  values.  By 
the  discussion  in  the  previous  section,  the  Nyquist  plot  of  an  optimal  compensator  will  therefore 
lie  on  each  of  these  shifted  boundaries  at  the  corresponding  w-value.  In  particular,  if  each  of 
the  shifted  boundaries  passes  through  the  point  (0,0),  then  no  compensation  is  necessary,  and 
an  optimal  compensator  has  been  found. 

After  viewing  this  display,  the  designer  is  given  the  option  to  do  nothing  or  to  design  a 
compensator  of  one  of  the  six  standard  forms  listed  in  the  previous  section.  If  the  design  requires 
a  pure  integrator  in  order  to  satisfy  specifications  on  steady-state  error,  then  the  designer  chooses 
option  6,  and  is  prompted  to  enter  the  value  of  K  in  dB  on  the  diagram  by  clicking  the  left 
mouse  button  at  any  point  on  the  K  dB  line. 

The  option  of  using  the  first  five  standard  forms  is  the  heart  of  OptComp,  and  this  is 
where  the  designer’s  skill  and  judgement  are  crucial.  However,  since  the  program  runs  quickly 
and  gives  clear  graphical  and  tabular  output,  poor  choices  are  immediately  obvious  and  easy  to 
discard.  Thus  experience  and  confidence  are  rapidly  acquired. 

After  viewing  the  boundaries,  the  designer  must  decide  on  the  lowest  order  standard  com¬ 
pensator  whose  Nyquist  plot  seems  most  likely  to  pass  close  to  them  (in  the  correct  order  of 
frequencies).  He  or  she  is  then  prompted  to  use  the  mouse  to  enter  on  the  diagram  (in  order) 
n  points  through  which  the  Nyquist  plot  of  the  compensator  should  pass  at  the  n  values  of  w. 
For  optimal  compensation,  these  points  should  be  on  the  respective  boundaries,  and  the  judge¬ 
ment  comes  firstly  in  choosing  the  form  of  compensator,  and  secondly  in  selecting  realizable 
points  close  to  the  boundaries.  The  program  fits  the  best  compensator  of  the  desired  form,  and 
re-displays  the  previous  diagram  with  the  addition  of  the  chosen  points  (green  circles)  and  the 
Nyquist  plot  of  the  fitted  compensator  (red  curve  and  stars).  The  green  circles  and  red  stars 
are  all  labelled  with  the  appropriate  w-values.  The  designer  can  thus  assess  visually  how  good 
his  or  her  choices  were,  and  how  effective  the  compensation  is.  (A  warning  message  appears 
if  the  fitted  compensator  is  not  minimum  phase.)  More  precise  information  is  then  given  in 
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tabular  form,  displaying  the  minimum  compensation  required  at  each  w-value,  firstly  for  the 
uncompensated  plant  and  then  with  the  new  provisional  compensator. 

It  is  now  easy  for  the  designer  to  decide  whether  the  new  compensator  is  satisfactory.  If 
not,  then  he  or  she  can  either  quit,  or  repeat  the  step  of  choosing  a  form  of  compensator  and 
selecting  the  points  on  the  diagram  through  which  its  Nyquist  plot  should  pass.  To  facilitate  this 
choice  at  second  and  later  attempts,  the  diagram  also  includes  the  red  stars  actually  achieved 
by  the  previous  unsatisfactory  attempt.  In  this  way,  the  designer  is  guided  into  a  better  choice 
the  next  time.  With  a  little  practice,  one  quickly  becomes  familiar  with  the  standard  Nyquist 
plots,  and  learns  what  choices  are  realizable  in  practice. 

If  the  designer  is  satisfied  with  the  compensation  achieved  by  this  form  of  compensator, 
then  the  gain  and  the  coefficients  of  numerator  and  denominator  are  displayed  and  written  to 
to  the  diary.  The  coefficients  appear  as  vectors,  so,  for  example,  the  form  4  compensator 

123(s^  +  4s  +  5) 

-b  6s  +  10 

would  appear  as  having  gain  123,  numerator  14  5,  and  denominator  1  6  10.  The  previously 
displayed  table  of  minimum  errors,  before  and  after  this  compensation  step,  is  also  written  to 
the  diary.  Next,  the  Nyquist  plot  of  the  plant  and  composite  boundaries,  with  which  the  process 
starts,  is  displayed  again,  now  also  including  the  Nyquist  plot  of  the  compensated  plant  (in 
green,  with  labelled  circles  at  the  w-values).  A  warning  message  appears  if  the  compensated 
plant  is  likely  to  have  closed  loop  instability  (i.e.,  if  its  Nyquist  plot  crosses  the  ±180°  phase 
line  above  the  0  dB  line). 

The  designer  can  then  quit  (if  this  compensator  has  adequately  met  the  design  specifica¬ 
tions)  or  repeat  the  process,  starting  with  the  already  compensated  plant,  to  increase  the  order 
of  the  compensator,  i.e.,  to  include  an  additional  compensator  of  one  of  the  standard  forms. 

5.  Comments  and  conclusions 

We  believe  that  the  procedures  in  OptComp  can  give  the  designer  all  the  transparency  of 
QFT  and  simplify  the  art  of  loop-shaping,  while  freeing  him  or  her  from  laborious  calculation. 
Experienced  designers,  familiar  with  traditional  loop-shaping  techniques,  may  need  a  period  to 
adjust  to  the  new  approach  of  working  directly  with  the  compensator  loop.  It  might  be  valuable 
to  try  the  program  with  students,  who  do  not  have  prejudices  about  other  methods,  to  see 
how  quickly  they  can  acquire  skill  in  loop-shaping  using  OptComp.  Experimentation  with  the 
program  will  also  have  the  advantage  of  suggesting  where  it  can  be  improved,  because  many  of 
its  features  have  come  from  our  own  trials  on  sample  data. 

A  full  listing  of  OptComp. m  appears  in  Appendix  3,  and  includes  many  comments  to  clarify 
the  structure  of  the  program  and  its  various  procedures.  The  curve-fitting  step  is  certainly  not 
the  best  possible,  although  the  speedy  trial-and-error  method  allows  points  giving  a  poor  fit  to 
be  simply  discarded  and  replaced  by  better  choices.  At  this  preliminary  stage,  the  compensator 
coefficients  are  found  by  linearizing  the  equation  in  two  ways,  using  Matlab’s  least  squares 
procedure  each  time,  then  taking  the  average.  Use  of  a  sophisticated  non-linear  optimization 
routine  might  give  closer  fits. 
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Appendix  1.  Glossary 


compensator  A  rational  function  of  s,  which  has  real  coefficients,  and  in  which  the 

degree  of  the  numerator  is  not  greater  than  that  of  the  denominator. 
A  transfer  function  is  multiplied  by  a  compensator  to  make  it  meet 
predetermined  specifications. 

controller  Alternative  term  for  compensator. 

m-file  An  executable  Matlab  file,  with  extension  .m. 

minimum  phase  A  minimum  phase  compensator  is  one  with  no  poles  or  zeros  in  the  right 

half-plane. 

Nichols  chart  Representation  of  complex  numbers  in  terms  of  their  phase  angles  (in 

degrees)  and  gain  magnitudes  (in  dB,  i.e.,  on  a  logarithmic  scale).  In  the 
programs,  each  point  in  the  Nichols  chart  is  treated  as  a  complex  number 
phase  -|-j  gain.  Except  for  the  scales  on  the  axes,  the  Nichols  chart 
representation  of  a  complex  number  is  equal  to  j  ln(z),  so  multiplication 
in  the  complex  plane  corresponds  to  addition  in  the  Nichols  chart.  The 
m-files  nic.m  and  denic.m  respectively  transform  complex  numbers  to 
and  from  their  Nichols  chart  representation,  using  the  formulae 
nic(z)  =  ^  arg(z)  +  20ilogio  \z\  and 
denic(tu)  =  exp(jY|5' Re(ti>)  -f  0.05  In  lOlm(tn)). 

Nyquist  plot  The  set  of  points  obtained  by  transforming  the  imaginary  axis  by  a  given 

function.  The  Nyquist  plot  of  /  is  the  set  {/(jw)  |  — oo  <  w  <  w}.  If 
f{z)  —  f{z),  then  the  plot  is  symmetrical,  and  only  the  portion  for  a;  >  0 
is  usually  plotted.  For  our  purposes,  all  Nyquist  plots  are  of  this  form 
and  are  in  the  Nichols  chart,  i.e.  {nic(/(ju>))  |  w  >  0}. 

order  For  a  rational  function,  the  sum  of  the  degrees  of  numerator  and  de¬ 

nominator. 

pole  excess  For  a  rational  function,  the  difference  between  the  degrees  of  denomina¬ 

tor  and  numerator. 

QFT  Quantitative  Feedback  Theory. 
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Appendix  2.  A  typical  OptComp  session 


The  following  is  the  complete  transcript  of  a  session  with  OptComp,  using  the  data  loaded 
from  OptCompl.  The  total  elapsed  time  is  19  minutes,  most  of  which  was  spent  waiting  for  the 
screen  dumps  of  the  graphical  displays. 

»  optcompl 

This  pleoit  can  be  compensated  with  a  linear/linear  compensator. 

»  opt comp 

Remember  omega  plant  bndry  must  be  defined  beforehand. 

Important  output  is  recorded  in  a  diary. 

Please  enter  0  to  quit  or  1  to  continue:  1 
12-Apr-96 
Starting  time: 

14  26 

The  next  diagram  gives  the  original  plant  values 

and  the  composite  boundaries  at  the  values  of  omega  indicated. 

(Press  einy  key.) 


Minimum  compensation  required  for  optimal  trainsfer  function 


omega 

dB 

angle 

0.4000 

5.2800 

-29.3600 

1.0000 

11.6200 

-13.6800 

3.0000 

16.0400 

37.9500 

10.0000 

20.1700 

90.0000 

40.0000 

29.5700 

-5.2100 

Optimal  compensator  values  must  lie  on  each  of  the  curves  shown  in 
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Please  enter  0  to  quit, 

1  for  constant/linear  (1st  order  lag) 

2  for  linear/linear  (1st  order  lead-lag) 

3  for  linear /quadratic 

4  for  quadratic/quadratic  (2nd  order  lead-lag) 

5  for  constant  over  quadratic  (2nd  order  lag) 

6  for  pure  integrator  (K/s) 


Use  the  mouse  to  choose  n  points  near  the  n  boundaries  shown 
(Press  any  key.) 

40, - 


Minimum  compensation  required  for  optimal  transfer  function 
Previously  With  this  compensator 


omega 

dB 

angle 

dB 

angle 

0.4000 

5.2800 

-29.3600 

-2.8379 

15.4838 

1.0000 

11.6200 

-13.6800 

-4.7049 

4.4256 

3.0000 

16.0400 

37.9500 

-3.0035 

-26.8701 

10.0000 

20.1700 

90.0000 

-0.1302 

12.8050 

40.0000 

29.5700 

-5.2100 

2.6032 

-3.2345 

Please  enter  0  if  satisfied  or  1  to  replace  this  compensator  or  quit:  1 

Optimal  compensator  values  must  lie  on  each  of  the  curves  shown  in 
the  next  diagram  at  the  values  of  omega  indicated.  (Press  ainy  key.) 


Please  enter  0  to  quit, 

1  for  constant/linear  (1st  order  lag) 

2  for  linear/linear  (1st  order  lead-lag) 

3  for  linear /quadratic 

4  for  quadratic/quadratic  (2nd  order  lead-lag) 

5  for  constant  over  quadratic  (2nd  order  lag) 

6  for  pure  integrator  (K/s) 

2 

Use  the  mouse  to  choose  n  points  near  the  n  boundaries  shown. 
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(Press  any  key.) 


Minimum  compensation  required  for  optimal  transfer  function 
Previously  With  this  compensator 


omega 

dB 

angle 

dB 

angle 

0.4000 

5.2800 

-29.3600 

-0.4857 

9.0292 

1.0000 

11.6200 

-13.6800 

0.8641 

5.5949 

3.0000 

16.0400 

37.9500 

-3.0216 

2.5373 

10.0000 

20.1700 

90.0000 

0.1979 

-10.5947 

40.0000 

29.5700 

-5.2100 

1.0117 

-4.1553 

Please  enter  0  if  satisfied  or  1  to  replace  this  compensator  or  quit:  0 
Compensator  gain 
27.9524 

Compensator  numerator 
1.0000  1.0976 

Compensator  denominator 
1.0000  9.1885 

Minimum  compensation  required  for  optimal  trainsfer  function 
Previously  With  this  compensator 


omega 

dB 

cuigle 

dB 

angle 

0.4000 

5.2800 

-29.3600 

-0.4857 

9.0292 

1.0000 

11.6200 

-13.6800 

0.8641 

5.5949 

3.0000 

16.0400 

37.9500 

-3.0216 

2.5373 

10.0000 

20.1700 

90.0000 

0.1979 

-10.5947 

40.0000 

29.5700 

-5.2100 

1.0117 

-4.1553 
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Appendix  3.  Files  on  disk 

The  enclosed  disk  (in  1.44  MB  3|  inch  DOS  format)  contains  six  m-files,  the  ifiain  program 
(OptComp.m),  two  small  function  files  (nic.m  and  denic.m)  for  conversion  to  and  from  the 
Nichols  chart,  and  three  sample  datafiles  (OptCompl.m,  0ptComp2.m,  and  0ptComp3.m),  any 
one  of  which  can  be  run  to  provide  input  data  for  OptComp.  A  listing  of  the  file  OptComp.m 
follows. 

•/.Program  OPTCOMP.M 

'/.This  is  an  interactive  program  for  designing  a  compensator  (or  controller) 

'/.  for  a  given  plant. 

'/.Input  data  (to  be  defined  beforehand)  are: 

'/.  omega  (1  by  n  vector  of  emgular  velocities) 

'/.  plant  (1  by  n  vector  of  nominal  plant  values  in  Nichols  chart, 

'/.  one  value  for  each  value  of  omega) 

'/.  bndry  (m  by  n  matrix  of  composite  boundary  points, 

'/.  m  boundary  points  for  each  value  of  omega) 


'/. - ; - PREAMBLE - 

disp(’ Remember  omega  plant  bndry  must  be  defined  beforehand.’) 

disp(’ Important  output  is  recorded  in  a  diary.’) 

itrt  =  input ( ’Please  enter  0  to  quit  or  1  to  continue:  ’); 


if  itrt  <  1 
return 
end 

•/ - INTRODUCTION 

diary  on 
disp(date) 

temp  =  round ( clock)  ,* 
temp  =  temp([0  00  1  10]); 
disp( ’Starting  time:’) 
disp(temp) 
diary  off 

[m,n]  =  size(bndry); 
for  1  =  2:n 


if  abs  (real  (plant  (1))  -  real  (plant  (1-1)))  >  200.0  '/.'/.This  extends  the 
if  real  (plant  (1))  >  real  (plant  (1-1))  '/.Nichols  chart,  if 

plant (1)  =  plant (1)  -  360.0;  '/.necessary,  to  make 

else  '/.the  plant  locus 

plant  (1)  =  plant  (1)  +  360.0;  '/.continuous 

end  */, 

end  •/, 

end  '/.I  n 
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J  =  real (plant)  >  180.0; 
temp  =  sum( J) ; 
if  temp  >  0.5*n 
plant  =  plant  -  360.0; 
else 

J2  =  real (plant)  <  -180.0; 
temp2  =  sum( J) ; 
if  temp2  >  0.5*n 
plant  =  plant  +  360.0; 
end 
end 

for  1  =  l:n 
for  k  =  2:m 


*/,*/, This  ensures  that 
'/the  majority 
*/,of  points  lie  between 
'/the  -180  and  +180 
'/degree  phase  lines 
'/ 

'/ 

'/ 

'/ 


if  abs(real(bndry(k,l))  -  real(bndry(k-l ,1)))  >  200.0  '/This  extends 

if  real(bndry(k,l))  >  real (bndry(k- 1,1))  '/the  Nichols  chart, 

bndry(k,l)  =  bndry(k,l)  -  360.0;  '/if  necessary,  to 

else  '/make  the 


bndry(k,l)  =  bndry(k,l)  +  360.0; 
end 


'/boundaries 

'/continuous 


end 

end  '/  k 

J  =  real(bndry(: ,1))  >  180.0; 
temp  =  sum(J) ; 
if  temp  >  0.5*m 

bndry(:,l)  =  bndry(:,l)  -  360.0; 
else 

J2  =  real(bndry( : ,1))  <  -180.0; 
temp2  =  sum(J) ; 
if  temp2  >  0.5*m 
bndry(:,l)  =  bndry(:,l)  +  360.0; 
end 


'/ 

'/ 

'/'/This  ensures  that 
'/the  majority 
'/of  points  lie  between 
'/the  -180  and  +180 
'/degree  phase  lines 
'/ 

'/ 

'/ 

'/ 

'/ 


end  •/•/ 

end  '/I  •/ 

disp(’The  next  diagram  gives  the  original  plant  values’) 
disp(’and  the  composite  boundaries  at  the  values  of  omega  indicated.’) 
disp(’ (Press  any  key.)’) 


pause 

plot (real (bndry) ,  imag(bndry) , ’-w’ , . . . 

real (plant ) ,imag (plant) , ’-b’ , real (plant ) ,imag (plant) , ’+b’) 
grid 

xlabeK ’Phase  (degrees)’) 
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ylabelC’Gain  (dB)’) 
for  1  =  l:n 

text (real (bndry (1,1)), imag (bndry (1,1)), num2str ( omega ( 1 ) ) ) 
t  ext (real (bndry (m , 1 ) ) , imag (bndry (m , 1 ) ) , num2st r ( omega ( 1 ) ) ) 
text (real (plant (1) ) , imag (plant (1) ) ,num2str (omega(l) ) ) 
end 
pause 

plant  1  =  plaint; 
ommin  =  0.5*min(omega) ; 
ommax  =  2.0*max(omega) ; 

omg  =  exp(linspace(log(ommin) ,log(ommax) ,100)) ; 


THE  MAIN  LOOP,  DESIGNING  ONE  STANDARD  COMPENSATOR 


while  itrt  >  eps  ‘/this  loop  goes  to  the  end  of  the  program  without  indentation 
plntmtrx  =  ones(m,l)*plantl; 

diff  =  nic ( deni c (bndry  -  plntmtrx));  The  difference 

for  1  =  l:n  '/.between  boundaries 

for  k  =  2:m  '/.and  plant  values, 

if  abs(real(diff (k,l))  -  real(diff (k-1,1)))  >  200.0'/,i.e.  curves  of  optimal 
if  real  (diff  (k,l))  >  real  (diff  (k-1,1))  '/.compensation  values 

diff(k,l)  =  diff(k,l)  -  360.0;  '/. 


else 

diff(k,l)  =  diff(k,l)  +  360.0; 
end 
end 

end  '/.  k 

J  =  reaKdiff  (:  ,1))  >  180.0; 
temp  =  sum( J) ; 
if  temp  >  0.5*m 
diff(:,l)  =  diff(:,l)  -  360.0; 
else 

J2  =  reaKdiff  (:  ,1))  <  -180.0; 
temp2  =  sum(J) ; 
if  temp2  >  0.5*m 

diff(:,l)  =  diff(:,l)  +  360.0; 
end 
end 

end  '/.I 

rediff  =  reaKdiff ); 
imdiff  =  imag (diff); 
temp2  =  abs(denic(diff )-l .0) ; 


'/.(Made  continuous 
'/.  by  extending  beyond 
'/.  +/-180  degrees 
'/.  if  necessary.) 

•/.•/. 

'/.'/.This  ensures  that 
'/.the  majority 
'/.of  points  lie  between 
'/.the  -180  and  +180 
'/.degree  phase  lines 
'/. 

•/. 

y. 

y. 

y. 

y.y. 


'/.'/.  This  finds  one 
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temp3  =  ones(m,  l)*min(temp2) ;  '/,  point  on  each  curve 

J  =  abs(temp2-temp3)<eps  ;  '/,  where  the  required 

J2  =  cumsuin(J)<2;  '/,  compensation 

J  =  J.*J2;  y,  has  smallest  magnitude, 

error  =  nic(denic(diff  (J))) ; 
diary  on 

disp( ’Minimum  compensation  required  for  optimal  transfer  function’) 

disp(’  omega  dB  angle’) 

disp([omega  ;  imag(error) . ’  ;  real(error) . ’] . ’) 

diary  off 

rpt  =  2; 

y. - - INNER  LOOP,  ALLOWING  REPEATED  ATTEMPTS - 

while  rpt>eps 

disp( ’Optimal  compensator  values  must  lie  on  each  of  the  curves  shown  in’) 
dispC’the  next  diagreim  at  the  values  of  omega  indicated.  (Press  any  key.)’) 
pause 

plotCrediff,  imdiff,’-w’) 
grid 

xlabeK ’Phase  (degrees)’) 
ylabeK’Gain  (dB)’) 
for  1  *  l:n 

text(rediff (1,1) .imdiff (1,1) ,num2str(omega(l))) 
text(rediff (m,l) , imdiff (m,l) ,num2str(omega(l))) 
end 

if  rpt  <  2 
hold 

plot (real (Gatt) ,imag(Gatt) , ’*r’) 
for  1  =  l:n 

text (real (Gatt (1)) ,imag(Gatt(l)) ,num2str(omega(l))) 
end 
hold 


end 

pause 

disp(’ Please  enter  0 


disp(’  1 
disp(’  2 
disp(’  3 
disp(’  4 
disp(’  5 
disp(’  6 


nc  =  input ( ’  ’ ) ; 


to  quit , ’ ) 

for  constant/linear  (1st  order  lag)’) 
for  linear/linear  (1st  order  lead-lag)’) 
for  linear/quadratic’) 

for  quadratic/quadratic  (2nd  order  lead-lag)’) 
for  constant  over  quadratic  (2nd  order  lag) ’ ) 
for  pure  integrator  (K/s)’) 
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nc  =  roimd(nc) ; 
if  nc  <  eps 
rpt  =  0; 
else 

if  nc  ==  6 

dispC’The  phase  is  constamt  at  -90  degrees.’) 

dispC’Use  the  mouse  to  choose  the  numerator  (i.e.  the  gain  at  omega  =  1)’) 

disp(’ (Press  any  key.)’) 

pause 

[temp,  K]  =  ginput(l); 

Gnum  =  [K]  ; 

Gden  =  [10]; 
else 

dispC’Use  the  mouse  to  choose  n  points  near  the  n  boundaries  shown.’) 

disp(’ (Press  auiy  key.)’) 

pause 

[rGreq,iGreq]  =  ginput(n); 

Greq  =  rGreq  +  i*iGreq; 
s  =  i* (omega’ ); 
sinv  =  ones(n,l) ./s; 
rhs  *  denic(Greq); 
if  nc  =®1 


'/This  fits  a  constant  over  line> 
Pcomp  =  [ones(n,l)  -rhs]; 

Qcomp  =  rhs.*s; 

P  =  [real (Pcomp) ; imag(Pcomp)]  ; 

Q  =  [real (Qcomp) ;imag(Qcomp)] ; 
cffsl  =  P\Q; 

Pcomp  =  [sinv  -sinv.*rhs]; 
Qcomp  =  rhs; 

P  =  [real (Pcomp) ; imag(Pcomp)] ; 

Q  =  [real (Qcomp) ;imag(Qcomp)]  ; 
cffs2  =  P\Q; 

cffs  =  0.5*(cffsl  +  cff s2) ; 

Gnum  =  [cff s(l)] ; 

Gden  =  [1  cffs (2)] ; 
end  '/,  of  case  nc  =  1 
if  nc  ==  2 

'/This  fits  a  linear  over  linear 
Pcomp  =  [s  ones(n,l)  -rhs]; 
Qcomp  =  rhs.*s; 


compensator  '/'/The  curve-fitting 

'/routines  are  fairly 
'/naive  - 

'/linearize  the  problem 
'/by  cross-multiplying, 
'/then  break  into  real 
'/eind  imaginary  parts, 
'/and  use  Matlab’s 
'/least  squares  fit. 
'/Then  divide  through 
'/the  highest  power  of  s 
'/and  repeat . 

'/Finally,  average  the 
'/two  answers . 

'/'/ 

compensator 
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P  =  [real(Pcomp) ;imag(Pcomp)]  ; 

Q  =  [real(Qcomp) ;imag(Qcomp)] ; 
cffsl  =  P\Q; 

Pcomp  =  [ones(n,l)  sinv  -sinv.*rhs]; 

Qcomp  =  rhs; 

P  =  [real (Pcomp) ; imagCPcomp)] ; 

Q  =  [real (Qcomp) ; imag(Qcomp)] ; 
cffs2  =  P\Q; 

cffs  =  0.5*(cffsl  +  cff s2) ; 

Gnm  =  [cffs(l)  cffs(2)]; 

Gden  =  [1  cffs (3)]; 
end  */,  of  case  nc  =  2 
if  nc  ==  3 

'/.This  fits  a  linear  over  quadratic  compensator 
Pcomp  =  [s  ones(n,l)  -rhs.*s  -rhs]; 

Qcomp  =  rhs.*s.“2; 

P  =  [real (Pcomp) ;imag(Pcomp)] ; 

Q  =  [real (Qcomp) ;imag(Qcomp)] ; 
cffsl  *  P\Q; 

Pcomp  =  [sinv  sinv. “2  -sinv.*rhs  -rhs.*sinv.*2] ; 

Qcomp  *  rhs; 

P  *  [real (Pcomp) ;imag(Pcomp)] ; 

Q  =  [real (Qcomp) ;imag(Qcomp)] ; 
cffs2  =  P\Q; 

cffs  =  0.5* (cff si  +  cff s2) ; 

Gnum  =  [cffs(l)  cffs(2)]; 

Gden  =  [1  cffs (3)  cffs (4)]; 
end  '/,  of  case  nc  =  3 
if  nc  ==  4 

'/.This  fits  a  quadratic  over  quadratic  compensator 
Pcomp  =  [s.“2  s  ones(n,l)  -rhs.*s  -rhs]; 

Qcomp  =  rhs.*s.~2; 

P  =  [real (Pcomp) ;imag(Pcomp)] ; 

Q  =  [real (Qcomp) ;imag(Qcomp)] ; 
cffsl  =  P\Q; 

Pcomp  =  [ones(n,l)  sinv  sinv. “2  -sinv.*rhs  -rhs.*sinv.“2] ; 
Qcomp  =  rhs; 

P  =  [real (Pcomp) ;imag(Pcomp)] ; 

Q  =  [real (Qcomp) ;imag(Qcomp)] ; 
cffs2  =  P\Q; 

cffs  =  0.5*(cffsl  +  cf f s2) ; 


19 


Gnum  =  [cffs(l)  cffs(2)  cffs(3)]; 

Gden  =  [1  cffs(4)  cffs(5)]; 
end  '/,  of  case  nc  ==  4 
if  nc  ==  5 

'/This  fits  a  consteint  over  quadratic  compensator 
Pcomp  =  [ones(n,l)  -rhs.*s  -rhs] ; 

Qcomp  =  rhs.*s.“2; 

P  =  [real (Pcomp) ;imag (Pcomp)] ; 

Q  =  [real (Qcomp) jimag (Qcomp) ] ; 
cffsl  =  P\Q; 

Pcomp  =  [sinv.“2  -sinv.*rhs  -rhs.*sinv."2] ; 

Qcomp  =  rhs; 

P  =  [real (Pcomp) ;imag(Pcomp)] ; 

Q  =  [real (Qcomp) ; imag(Qcomp)] ; 
cffs2  =  P\Q; 

cffs  =  0.5* (cffsl  +  cff s2) ; 

Gnum  =  [cffs(l)]; 

Gden  =  [1  cffs(2)  cffs(3)]; 
end  y,  of  case  nc  ==  5 
end  y,  of  case  nc  >  0 
if  Gnum  >=  -eps 
else 

disp( ‘ *****************WARNING**************************** ’ ) 
disp(’This  compensator  has  a  zero  in  the  right  half -plane .’ ) 
disp( ’ *****************WARNING**************************** ’ ) 
end 

if  Gden  >=  -eps 
else 

disp ( ’ *****************WARNING**************************** ’ ) 
disp(’This  compensator  has  a  pole  in  the  right  half -plane .’ ) 
disp(’  ************=t!****WARNING****************************  ’ ) 
end 

Gatt  =  nic (polyval (Gnum, i*omega) ./polyval (Gden, i*omega)) ; 

Gcrv  =  nic (polyval (Gnum, i*omg) ./polyval (Gden, i*omg)) ; 
if  nc  "  6 
Greq  =  Gatt; 
end 

plot(rediff ,  imdiff , ’-w’ , real (Greq) ,imag(Greq) , ’og' , real (Gatt) , . 
imag(Gatt) , ’*r’ ,real(Gcrv) ,imag(Gcrv) , ’-r’ ) 

grid 

xlabeK ’Phase  (degrees)’) 
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ylabeK’Gain  (dB)’) 
for  1  =  l:n 

text(rediff (1,1) .imdiff (1,1) ,nvim2str(omega(l))) 
text(rediff (m,l) , imdiff (m,l) ,num2str(omega(l))) 
text (real (Greq(l)) ,imag(Greq(l)) ,num2str(omega(l))) 
text (real (Gatt (1) ) , imag(Gatt (1) ) ,num2str (omega(l) ) ) 
end 
pause 

newdiff  =  diff  -  ones(m,l)*Gatt;  This  computes  how  much 
temp2  =  abs(denic(newdiff  )-l .0) ;  */,  extra  compensation  is 

temp3  =  ones(m,l)*min(temp2) ;  */,  needed,  if  the  current 

J  =  abs(temp2-temp3)<eps  ;  '/,  compensator  is  used. 

J2  =  cumsum(J)<2;  '/,  (cf  diff  and  error  above) 

J  =  J.*J2;  •/, 

newerr  =  nic(denic(newdiff  (J))) ; 

disp( ’Minimum  compensation  required  for  optimal  transfer  function’) 
disp(’  Previously  With  this  compensator’) 

disp(’  omega  dB  angle  dB  auigle’) 

disp( [omega’  imag( error)  real (error)  imag (newerr)  real (newerr) ] ) 
rpt  =  input (’Please  enter  0  if  satisfied  or  1  to  replace  this  compensator  or 
quit:  ’); 

end  */,  of  else  nc  >=  eps 
end  '/,  of  while  rpt  >  eps 


I - end  of  inner  loop - - - 

'/. - FINAL  SECTION  OF  THE  MAIN  LOOP 


if  nc  >  eps 
diary  on 

if  abs(Gnum(l))  >  eps 
K  =  Gnum(l) ; 

Gnum  =  Gnum/K; 

disp(’ Compensator  gain’); 

disp(K) 


*/,This  section  applies  only  if 
'/,a  standard  compensator  has 
'/been  designed. 

'/It  gives  the  coefficients, 
'/plots  the  boundaries  eind 
'/compensated  and  uncompensated 
'/plant  loci, 

'/and  tabulates  the  minimum 
'/compensation  necessary,  before 
'/and  after  introducing  this 
'/compensator. 


disp( ’Compensator  nvimerator’)  '/compensation  necessary,  before 

disp(Gnum)  '/and  after  introducing  this 

disp(’ Compensator  denominator’)  '/compensator. 
disp(Gden) 

disp( ’Minimum  compensation  required  for  optimal  trauisfer  function’) 
disp(’  Previously  With  this  compensator’) 
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disp(’  omega  dB  angle  dB  eoigle’) 

disp( [omega’  imag(error)  real(error)  imag(newerr)  real(newerr)] ) 
diary  off 

plantl  =  nic(denic(plantl  +  Gatt)); 
temp  =  0; 
for  1  =  2:n 

if  abs (real (plantl (1))  -  real (plantl (1-1)))  >  200.0 
if  temp  <  1 

if  imag(plantl(l))  >  0 

disp( ’ *****************WARNING**************************’ ) 
disp(’The  closed  loop  transfer  function  may  be  unstable.’) 
disp( ’*****************WARNING**************************’ ) 
temp  =  1; 
end 
end 

if  real (plantl (1))  >  real (plantl (1-1)) 
plantl(l)  =  plauitld)  -  360.0; 
else 

plantl(l)  =  plantl(l)  +  360.0; 
end 
end 
end  ‘/,1 


J  =  real(plantl)  >  180.0; 
temp  =  sum(J) ; 
if  temp  >  0.5*n 

plantl  =  plantl  -  360.0; 
else 

J2  =  real(plantl)  <  -180.0; 
temp2  =  sum(J); 
if  temp2  >  0.5*n 

plantl  =  plantl  +  360.0; 
end 


*/,7,This  ensures  that 
*/,the  majority 
'/,of  points  lie  between 
'/.the  -180  eind  +180 
'/.degree  phase  lines 
'/. 

'/. 

'/. 

'/. 


end  1% 

disp(’The  next  diagram  gives  the  original  and  compensated  plant  values’) 
disp(’eind  the  composite  boundaries  at  the  values  of  omega  indicated.’) 
disp(’ (Press  any  key.)’) 


pause 

plot (real (bndry) ,  imag(bndry) , ’-w’ .... 

real (plantl) ,imag(plantl) , ’-g’ .real (plant 1) ,imag(plantl) , ’og’ , . . . 
real(plant)  ,imag(plant) ,  ’-b’  .reaKpleint)  ,imag(plant) ,  ’+b’) 
grid 
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xlabelC ’Phase  (degrees)’) 
ylabeK’Gain  (dB)’) 
for  1  =  l:n 

t  ext (real (bndry (1,1)), imag (bndry (1,1)), num2str ( omega ( 1 ) ) ) 
text (real (bndry(m,l)) , imag (bndry (m,l)) ,num2str(omega(l))) 
text  (real  (plaint  1  (1) ) ,  imag(plantl  (1) )  ,num2str  (omega(l) ) ) 
t  ext  (  real  (plant  ( 1 ) ) ,  imag  (plaint  ( 1 ) ) ,  num2  str(omega(l))) 
end 
pause 

end  */,of  if  nc  >  eps 

y. - END  OF  FINAL  SECTION  OF  THE  MAIN  LOOP - 

disp( ’Please  enter  0  to  quit  or  1  to  design  a  higher  order  compensator:’) 

itrt  =  input ( ’  ’ ) ; 

end  '/,  of  while  itrt  >  eps 

'/. - - - END  OF  THE  MAIN  LOOP - 


diary  on 

temp  =  round(clock) ; 
temp  =  temp([0  00  1  10]); 
disp(’ Finishing  time:  ’) 
disp(temp) 
diary  off 
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